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[1] This study is devoted to the development of a semianalytical algorithm for the 
determination of the optical thickness, the liquid water path, and the effective size of 
droplets from spectral measurements of the intensity of light reflected from water clouds 
with large optical thickness. The probability of photon absorption by droplets in the visible 
and near-infrared spectral regions is low. This allows us to simplify and modify well- 
known asymptotic equations of the radiative transfer theory, taking into account the fact 
that the single scattering albedo is close to one. Modified asymptotic equations are used to 
develop the inverse algorithm. We also avoid the use of the Mie theory, applying 
parametrizations and geometrical optics results with account for wave corrections. The 
main advantage of the method proposed lays in the fact that the equations derived not only 
provide a valuable alternative to the numerical radiative transfer solution but also are much 
more simple than equations of a conventional asymptotic theory. This simplicity allows 


both the simplification of the cloud retrieval algorithm and, even more important, the 


insight into various factors involved in a retrieval scheme. 
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1. Introduction 


[2] Clouds play an important role in the Earth climate 
system [Kondratyev and Binenko, 1984; Liou, 1992]. The 
amount of radiation reflected by the Earth-atmosphere 
system into outer space does not only depend on the cloud 
cover and the total amount of condensed water in the Earth 
atmosphere. The size of droplets or crystals aef is also of 
importance. 

[3] The information about microphysical properties and 
spatial distributions of terrestrial clouds on a global scale 
can be achieved with satellite remote sensing systems only. 
Different spectrometers and radiometers [Bovensmann et 
al., 1999; Deschamps et al., 1994; King et al., 1992; 
Nakajima et al., 1998], deployed on space-based platforms, 
measure the angular and spectral distribution of intensity 
and polarization of the reflected solar light. Generally, the 
measured values depend on both geometrical and micro- 
physical characteristics of clouds. Thus, the inherent proper- 
ties of clouds can be retrieved (at least in principle) by 
solution of the inverse problem. The accuracy of the 
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retrieved values depends on the accuracy of measurements 
and the accuracy of the forward radiative transfer model. 

[4] In particular, it is often assumed that clouds can be 
represented by homogeneous and infinitely extended, in the 
horizontal direction plane-parallel slabs [Goloub et al., 2000; 
Han et al., 1994; King, 1981, 1987; Nakajima and King, 
1990; Nakajima et al., 1991; Rossow, 1989]. The range of 
applicability of such an assumption for real clouds is very 
limited as is readily shown by observations of light from the 
sky on the cloudy day. For example, the retrieved cloud 
optical thickness T apparently depends on the viewing 
geometry [Loeb and Davies, 1996; Loeb and Coakley, 
1998]. This, of course, would not be the case for an idealized 
plane-parallel cloud layer. However, the state-of-the-art radi- 
ative transfer theory and computer technology do not allow 
one to incorporate 3D effects into operational satellite 
retrieval schemes. As a result, cloud parameters retrieved 
should be currently considered as a rather coarse approxima- 
tion to reality. This justifies the development of new cloud 
retrieval algorithms, which can in principle have a lower 
accuracy than exact radiative transfer models, if one applies 
them to simulated radiance fields. However, they may yield 
very similar accuracies, when applied to real data. This is due 
to uncertainties of forward models for real atmospheres. 

[s] Note that even such limited tools produce valuable 
information on terrestrial clouds properties. For example, it 
was confirmed by satellite measurements that droplets in 
clouds over oceans are usually larger than over land [Han et 
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al., 1994]. This feature, for instance, is of importance for the 
simulation of the Earth’s climate [S/ingo, 1989]. 

[6] A new semianalytical retrieval algorithm for the cloud 
liquid water path and water droplets size determination has 
been developed and tested in this study. It is based on the 
asymptotical solution of the radiative transfer equation for a 
special case of disperse media with a large optical thickness. 
This solution was obtained by Germogenova [1963] for 
plane-parallel turbid slabs. Such an approach has already 
been used in a number of studies [Rozenberg et al., 1978; 
King, 1987]. The difference is that the asymptotical sol- 
utions are further simplified such that the inverse problem is 
reduced to the solution of a single transcendent equation. 
This allows us to speed up the retrieval process significantly 
with a minimal loss of accuracy of the retrieved parameters 
for simulated reflected light fields over optically thick 
clouds. We do not expect big differences between two 
models as far as real measurements are concerned. 

[7] The algorithm is restricted to the case of optically thick 
clouds. Thus, we assume that the optical thickness of clouds 
is larger than approximately 10. For small optical thick- 
nesses, the error increases. Thus, the conventional lookup 
table approach should be applied in the range of the optical 
thickness 0—10. In principle, if a high accuracy is not a 
primary target and one needs to have first order estimates, the 
algorithm can be also applied down to the optical thickness 5. 
So in this case the lookup table approach only needs to be 
applied to optically thin clouds with optical thickness t < 5. 
However, it should be stressed that the optical thickness is 
highly correlated with the geometrical thickness of clouds. 
For instance, it was shown [Feigelson, 1981] that clouds 
having T <5 are characterized by a geometrical thickness less 
than 200 m. It is difficult to expect that such thin clouds are 
horizontally homogeneous media. This requires taking the 
horizontal photon transport into account [Cahalan et al., 
1994; Platnick, 2001], which is not considered in standard 
retrieval procedures [Arking and Childs, 1985; Nakajima and 
King, 1990]. Due to this inhomogeneity, satellite cloud 
retrievals for optically thin clouds are troublesome even if 
the lookup table approach is used. Therefore, we will con- 
sider here only that part of our algorithm that is concerned 
with the case of optically thick clouds. 

[s] The paper is organized as follows. In section 2 we 
review the asymptotic theory of the radiative transfer [Rozen- 
berg, 1966; van de Hulst, 1980]. We also introduce here the 
main equation (33), which represents the reflection function 
of an optically thick weakly absorbing cloud layer over a 
underlying Lambertian surface with albedo A. The auxiliary 
functions in equation (33) are given in terms of elementary 
functions (except the reflection function of a semiinfinite 
nonabsorbing cloud (see, however, equation (5) for a special 
case of nadir measurements)). As a matter of fact, equation 
(33) is a very general one and can be used well outside of the 
framework of this work (e.g., for other weakly absorbing 
geophysical light scattering media like snow, ice, and foam). 

[9] Section 3 is devoted to the application of equation 
(33) to the semianalytical solution of the cloud inverse 
problem, namely for finding the effective radius of droplets 
and cloud liquid water path from a two-channel algorithm. 
Knowing this, the spectral cloud optical thickness and the 
number of cloud droplets in a unit column of a cloud 
volume can be easily obtained as well. Our solution of the 
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inverse problem is a semianalytical one only due to the 
necessity to solve a simple transcendent equation. This 
equation is solved using Brent’s method [Brent, 1973]. In 
practical terms, it does almost not influence the speed of the 
retrieval procedure compared to the case of a simple 
analytical function calculation. This feature could be of 
importance for operational algorithms, which deal with 
huge amounts of data. We furthermore study the measure- 
ment error propagation in the retrieval scheme (both in 
numerical and analytical forms) here. 

[10] In the conclusions we discuss advantages and dis- 
advantages of the modified asymptotic theory as applied to 
cloud parameter retrieval. Short appendices present simple 
analytical approximations for local optical characteristics of 
clouds like cloud phase function, absorption, and extinction 
coefficients. A simple and yet accurate parametrization for 
the cloud asymmetry parameter is also given there. Gen- 
erally, equations presented in this paper can be used in 
studies of other problems, well beyond the rather narrow 
topic of cloud satellite retrievals considered here. 


2. The Reflection Function 
2.1. The Visible Range 


[u] Our main aim in this section is to obtain a simple 
analytical result for the cloud reflection function in terms of 
elementary functions, which can be used for the semian- 
alytical satellite cloud retrieval algorithm development. 
Absorbing and nonabsorbing clouds are considered sepa- 
rately in sections 2.1 and 2.2, respectively. Section 2.3 is 
devoted to the derivation of new analytical equations for 
cloud spherical reflectance and total light absorptance in a 
cloud. Main factors involved in the retrieval procedure are 
discussed in section 2.4. 

[12] To start with, we note that the reflection function of a 
cloud R(ðo, 0, p) is defined as the ratio of the reflected light 
intensity Z(Ùo, 0, p) for the case of a cloud to that of an ideal 
Lambertian white reflector [Liou, 1992]: 


_ I(9o, Ù, p) 
Ro ea (1) 
where 
I* (o) = F cos ðo (2) 


is the intensity of light reflected from the ideally white 
Lambertian reflector, nF is the solar flux on the area 
perpendicular to the direction of incidence, Ùo is the solar 
angle, 0 is the observation angle and » is the relative 
azimuth between solar and observation directions. We use 
also: p = |cos I, po = cos Ùo. 

[13] It follows for the Lambertian ideally white reflector 
from equation (1): R = 1. This result does not depend on the 
viewing geometry by definition. Although clouds are white 
when looking from space, their reflection function R(Yo, Ù, 
y) is not equal to one. It depends on the viewing geometry. 
The results of calculations of the reflection function of an 
idealized semi-infinite nonabsorbing water cloud RO (Ùo, Ù, 
4) at the wavelength A = 650 nm and the nadir observation 
are presented in Figure 1. Calculations were performed for 
the gamma particle size distribution: 


f(a) = Dave, (3) 
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Figure 1. The reflection function of an idealized semiinfinite nonabsorbing cloud R20, Ùo, 0) at nadir 
observation obtained from the exact radiative transfer code [Mishchenko et al., 1999] and approximation 
(5) at wavelength A = 0.65 um and the effective radius of droplets a-y= 6 and 16 um. It is assumed that 
particles in a cloud are characterized by the gamma particle size distribution. 


s+1 

where 4 = TE (2) is the normalization constant and 
do is the mode radius. The half-width parameter s was equal 
to 6. This value is typical for terrestrial clouds. The effective 
radius, which is frequently used in cloud remote sensing 
studies, is defined as a ratio of the third to the second 
moment of the particle size distribution. In particular, for the 
distribution given by equation (3) we have: 


ae = a(t +3), (4) 


The values of the effective radius in Figure 1 were 6 and 
16 um. This covers the usual range of variability of the 
effective radius in natural water clouds. The function R(0,, 
Ù, p) can be smaller and larger than 1 depending on the 
incidence angle. This implies that, for particular viewing 
geometries, a cloud is even more reflective than the ideally 
white Lambertian surface. This is mostly due to peculia- 
rities of the phase function of cloudy media (e.g., in the 
backscattering (0 ~ Ùo, p ~ T) region) for comparatively 
thick clouds. It follows from Figure | that in the range of 
solar angles 30°—60° and nadir observation the reflection 
function of a water cloud is almost equal to the reflection 
function of an ideally white Lambertian reflector. It differs 
by not more than 10% from 1 for these geometries. 

[14] The phase functions p(8) (© is the scattering angle) of 
water clouds are depicted in Figure Al (see Appendix A). 


They are normalized by the following condition: 
5 fo p(8) sin6d0 = 1. Calculations were performed using 
the Mie theory at the wavelength of 0.65 um. The effective 
radii of water droplets were 4, 6, and 16 um. The value of s in 
equation (3) was set to be equal 6. The comparison of Figures 
1 and A1 indicates that the phase function is much more 
sensitive to the effective radius of droplets than the reflection 
function. We also see that all phase functions of water clouds 
pass a crossing point around the scattering angle 30°, having 
a value of approximately 2.27. The phase function is smaller 
than 1 at scattering angles larger than 43°. 

[15] The reflection function R? (ùo, Ù, p) can be repre- 
sented by the following simple approximate equation 
[Kokhanovsky, 2002a]: 


bı + by cos Ù cos Vo + p(8) 
4(cos 0 + cos Jo) 


Ri. (00,9, 9) = (5) 


where 6 = arccos (—cos Ù cos Ùo + sin 0 sin Ùo cos y) is the 
scattering angle, p(®) is the phase function of a cloudy 
medium (see Appendix A), bı and b2 are constants. 
Equation (5) obeys to the reciprocity principle [Zege et 
al., 1991]. It follows for nadir observations [Kokhanovshy, 
2002a, 2002b]: bı = 1.48, by = 7.76. The results of 
calculations of the function R% (Ùo, 0, p) with equation (5) 
are presented in Figure | along with data, obtained from 
exact radiative transfer computations. The comparison of 
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approximate and exact data shows that the accuracy of 
equation (5) is better than 2% at Vo < 85°. Constants bı and 
bə for other viewing geometries can be found using 
parametrizations of results obtained from exact radiative 
transfer codes [e.g., Mishchenko et al., 1999; Kokhanovsky 
and Rozanov, 2003]. 

[16] Equation (5) can be also used to find the reflection 
function of a finite cloud R(p, po, Y, T) by means of the 
following equation [Germogenova, 1963; van de Hulst, 
1980; Minin, 1988; Kokhanovsky, 2001]: 


R(p, Hos 9, T) = R% (u, Mo, P) — t(7)Ko(W)Ko(uo), (6) 


where 
1 


TE 
0.751(1—g)+a 


(7) 


is the global transmittance of a cloud, Ko(1) is the escape 
function. The escape function is defined via the solution of 
the characteristic integral equation [van de Hulst, 1980]. 
Note that parameters a and g in equation (7) are defined as 
follows [Sobolev, 1972]: 


1 


4 J Ko(u)p?dy, (8) 
0 
io ae 
g= 5 [ro sin 0 cos 6d0. (9) 
0 


[17] It should be pointed out that the escape function 
Ko() only weakly depends on the cloud microstructure and 
can be presented by the following simple equation [Sobolev, 
1972; Zege et al., 1991; Kokhanovsky, 2001]: 


(1 + 2p). (10) 


3 
Ko(u) =3 
The function Ko(\1), which was calculated with the exact 
radiative transfer code for g equal to 0.75, 0.85, and 0.9 in 
the case of Heney—Greenstein phase function, is presented 
in Figure 2. We see that Ko(\) does almost not depend on g 
at u > 0.2(0 < 78°). This is also the case for g = 0 and for 
the Mie-type phase functions. At the range of observation 
angles 0} = 80°—90° there is some dependence of the escape 
function on the microstructure of the cloud. However, the 
cloud top nonuniformity plays a role at such grazing 
observation angles. Therefore the problem cannot be solved 
in the plane-parallel layer approximation in this case 
anyway. 

[1s] The variability of Ko(\1) at u = 0.2—1.0 for different 
values of the average cosine of the scattering angle g = 
0.75—0.9 is well inside the 2% interval. This coincides 
with the error of equation (10) at u > 0.2, which is 
smaller than 2%. This confirms the wide range of applic- 
ability of equation (10) in cloud optics. Note that function 
(10) also describes the angular distribution of the diffused 
solar light transmitted by a thick cloud. In particular, for 
the ratio O of the diffused transmitted light intensity at the 
observation angle 60° (from the nadir) to that exactly at 
the nadir, we have from equation (10): © = 2/3, which 
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gives approximately 30% reduction. This can be easily 
checked experimentally. 

[19] The substitution of equation (10) into equation (8) 
yields: 


15 


=—# 1.07 
14 


(11) 


Q 


independently of the cloud microphysical parameters. It 
should be pointed out that the value of a, numerically 
calculated by King [1987], assuming the fair weather 
cumulus cloud model, is approximately 1.07 as well, i.e., in 
agreement with our estimation. King [1987] used the Mie 
theory to find the phase function of a cloudy medium. 
Yanovitskij [1997] found the same value of a for Heney— 
Greenstein phase functions with asymmetry parameters in 
the range 0.0—0.9. This supports the approximation using a 
fixed value of a in equation (7) given by (11), indepen- 
dently of cloud microphysical parameters. 

[20] The accuracy of equation (6) with account for (5), 
(7), (10), (11), which is used for approximate calculations of 
the reflection function of a finite cloud in visible, is 
illustrated in Figures 3a—3c. The error is less than 1% at 
T > 10 and A = 0.65 um at the nadir observation and a solar 
angle of 60° (see Figures 3a and 3b). This range of optical 
thicknesses coincides with the most frequently observed 
ones in water clouds, both using satellite and ground-based 
techniques [Trishchenko et al., 2001]. The small, almost 
constant error at T > 30 is mostly due to the error of 
approximation (5) for a semi-infinite cloud. Note that errors 
are negligibly small for all practical purposes for optically 
thick cloud fields. 

[21] The error depends on the solar angle (see Figure 3c). 
However, it does not exceed 5% at solar angles smaller than 
75° for values of tT > 10 and A = 0.65 um. The error 
increases for smaller cloud optical thicknesses (see Figures 
3a and 3b). It reaches 4% at T = 5 for the wavelength 0.65 
jum at the nadir observation and solar angle 60° (see Figure 
3b). It is highly sensitive to the solar angle at Ùo < 40° and T 
= 5 (see Figure 3c). The detailed study of the accuracy of 
equation (6) with account for equations (5), (7), (10), and 
(11) can be found elsewhere [Kokhanovsky and Rozanov, 
2003]. Kokhanovsky and Rozanov [2003] introduced cor- 
rections to the asymptotic result (6), which allowed to 
reduce the error to a level below 5% for all t > 5 for 
nongrazing solar zenith angles and nadir observation. How- 
ever, we will not consider this result here and assume the 
value of t > 10 for this study. 

[22] Equation (6) is readily modified to account for the 
light reflection from an Lambertian underlying surface 
[Sobolev, 1972]: 


A AT(W)T (Ho) 


R(t, bos P, T) = RCH, Ho, P, T) + (12) 
1 —Ar 


where Ê is the reflection function of a Lambertian surface- 
cloud system, A is the spherical albedo of an underlying 
Lambertian surface, which may depend on the wavelength, 
R= R(A=0), 


T(u) = Koln) (13) 


escape function 
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Figure 2. The escape function calculated with exact radiative transfer code for the Heney—Greenstein 
phase function at g = 0.75, 0.8, and 0.9 [ Yanovitskij, 1997] and with approximation given by equation (10). 
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Figure 3a. The dependence of the reflection function of a cloudy layer on the optical thickness 
according to equations (6) (at \ = 0.65 um) and (30) (at \ = 1.55 pm) at aef = 6 pm, Ù = 0°, 0 = 60° as 
compared to exact radiative transfer computations. 
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Figure 3b. The errors of approximations given by equations (6) (at A = 0.65 um) and (30) (at \ = 1.55 
jum) at aef = 6 pm, 0 = 0°, Do = 60° as compared to exact radiative transfer computations. 


is the diffuse transmittance of a cloud layer [Sobolev, 1972] andr neglected the direct solar light term in equation (12). This is 
is the spherical albedo of a cloud. Due to the energy conservation possible due to a large thickness of clouds under consideration. 
law we have: r = 1 — ż in visible, where we neglect possible [23] Finally, substituting equations (6) and (13) into 
small light absorption inside a cloud. Note that we have equation (12) we have as the reflection function of a 
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Figure 3c. The same as in Figure 3b, but as the function of the solar zenith angle at tT = 5 and T = 10. 
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Lambertian surface-cloud system: 


t(1 — A) 


T= AG — g&l). (14) 


Riu, Ho; P; T) = R% (p, Ho; $) > 


This formula with account for equations (5), (7), (10), and 
(11) will be used as a basis for the semianalytical cloud 
retrieval algorithm described in section 3. 


2.2. The Near-Infrared Range 


[24] Unfortunately, the relatively simple equation (6) 
cannot be applied to the calculation of the reflection 
function of a finite cloud in the near-infrared region of the 
electromagnetic spectrum. This is due to the presence of 
absorption bands of condensed water there. Alternatively, 
the following formula applies [Germogenova, 1963; van de 
Hulst, 1980; Sobolev, 1984; King, 1987]: 


mle?" 


R(M, Ho, P, T) = Roo (bs Ho, P) — To pe m K(W)K (Ho), (15) 


where y is the diffusion exponent, K(p) and Rẹ are the 
escape function and the reflection function of an absorbing 
semiinfinite medium with the same local optical character- 
istics as a finite layer under study, respectively. Equation 
(15) accounts for the influence of light absorption on the 
reflection function of clouds. Clearly, the reflection function 
decreases if additional absorbers are present in cloud 
droplets for the same phase function. Such a decrease is 
characteristic for the infrared channels, where water itself is 
a relatively strong absorber of solar radiation. 

[25] Constants / and m are defined by the following 
integrals [van de Hulst, 1980]: 


l=2 | K(p)i(~u)udp, (16) 
| 


m=2 f Pw nde, (17) 


i(n) being the angular distribution of light in deep layers of a 
cloud, where the so-called asymptotic regime takes place 
[Sobolev, 1972]. Functions i(), K(\), R? (p, Ho, 9) can be 
found from the solution of correspondent integral equations 
[Sobolev, 1972; van de Hulst, 1980]. 

[26] Functions Ralu, po, 9), K(p) and constants m and / 
have the following asymptotic forms when light absorp- 
tion by droplets is relatively small (single-scattering 
albedo wo — 1) [van de Hulst, 1980; Minin, 1988]: 


R(t hos 8) = Rolha, 8) = Ayar KOH Kolo), (18) 
x) =09( 1-20)" 


(19) 


(20) 
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(14) = (21) 
3(1 = g)’ 
and y 3(1 — wo) (1 — g) as the single scattering albedo 


wo — 1. Note that we have at wọ = 1: Rx, = RQ, K(w) = 
Ko(u), m = 0, l= 1. Taking correspondent limits, we obtain 
that equation (15) transforms into equation (6) in this case 
as it should be. Interestingly, as it follows from the 
asymptotic analysis [van de Hulst, 1980], the absorption 
reduces the cloud reflection function equally for all relative 
azimuths as wo — 1 (see equation (18)). 

[27] We underline that equations (18)—(21) follow from 
the asymptotic analysis of the radiative transfer equation. 
The integration of equation (18) with respect to all angles 
yields [Kokhanovsky, 2001]: r% = 1 — y, where y = 4,/3- 
and rœ is the spherical albedo of a weakly absorbing sem1- 
infinite cloud [van de Hulst, 1980]. This result is valid only as 
y — 0. Thus, the parameter y (as wo — 1 ) can be interpreted 
as a fraction of photons, absorbed in a weakly absorbing 
semi-infinite cloud for a diffuse illumination. It depends both 
on the single scattering albedo and the asymmetry parameter. 
Clouds having larger values of g, therefore, absorb more 
light. Larger values of g imply that photon scattering 
increases at small angles. Thus the photon path length before 
its escape from the medium is also increased. As a conse- 
quence, this results in increased light absorption in the 
medium. Media having different values of wo and g, but the 
same values of y, have the same values of r,, under the 
approximation considered here. The parameter y (divided by 
four and with gwo substituted by g) is called the similarity 
parameter [van de Hulst, 1980]. It is a useful parameter, 
describing the optical properties of clouds. Substituting y into 
equations (18)—(21) yields: 


Roo (Ut, Hos P) = RS (h, Hos 9) — ¥Ko(14) Ko (ho), (22) 
Ky) = Kow (1 - Z), (23) 
m = 2y, (24) 
and 
Pai (25) 


Equations (22)—(25) were derived assuming that wọ — 1. 
Alternatively, the right-hand sides of equations (22)—(25) 
give us the first terms of the expansion of correspondent 
functions with respect to y. The accuracy of equations 
decreases with the probability of photon absorption B = 1 — 
wo. The higher terms of the expansions (22)—(25) are not 
known or quite complex [Minin, 1988]. However, it has 
been shown that the following equations (see also equation 
(15)) approximately account for higher order terms [Rozen- 
berg, 1962; Zege et al., 1991; Kokhanovsky et al., 1998]: 


Ræ (p, Ho; p) = RY exp(—yu(p, Ho, 9)), (26) 


mK (u)K (po) = (1 — e™™)Ko(u)Ko (uo), (27) 
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l = exp(—ay), (28) 
where a viewing function is defined by 
Ko(u)Ko (ho) 
U({L, Lg, Y) = —— 29 


and does not depend on wo and 7 by definition. The viewing 
function only has a weak dependence on the microphysical 
parameters of clouds (e. g., the droplet size distribution). 
This fact follows from the low sensitivity of the functions 
R? (u, po, P), Kolu) in equation (29) to the microstructure of 
clouds (see Figures | and 2). 

[28] It should be stressed that the approximations (26)— 
(28), (5), (10), and (11) allow to avoid the solution of 
integral equations [van de Hulst, 1980; Yanovitskij, 1997] 
for the determination of functions R,,(|1, o, Y) and K(p) in 
equation (15). This is of a general importance. Note, that 
equations (26) and (28) transform to exact asymptotic 
results (22) and (25) as y — 0. 

[29] Substituting equations (27) and (28) into equation 
(15), we have: 

R(p, Hos P, T) = Rx (p, Ho; 9) = te ™ Ko()Ko(Ho); (30) 
where x = yt. Note that the value of exp(—x) describes the 
attenuation of light field in deep layers of a cloud [van de 
Hulst, 1980]. 

[30] It can be easily shown [Kokhanovsky, 2001] that the 
parameter ¢ in equation (30) gives the global transmittance 
for a nonabsorbing cloud field: 


sinh y 

~ sinh(ay + x) an) 
Importantly, equation (31) yields equation (7) at wo = 1. 

[31] Equation (30) was first proposed by Rozenberg 
[1962]. However, his derivation differs from that presented 
here. Also he proposed equation (31) with a = 1, which is 
not consistent with the exact asymptotic result, given by 
equation (7). 

[32] The range of applicability of equation (26) with 
respect to higher values of y can be readily extended using 
the following simple correction: 


Roo(Hs Hos P) = R3 (Hts Ho, P) exp(—Y(1 — ev)u(H, Ho, 9). (32) 
where c = 0.05 [Kokhanovsky, 2002b]. The value of c was 
obtained by the parametrization of exact radiative transfer 
calculations [Mishchenko et al., 1999]. 

[33] The accuracy of equations (30)—(32) with account for 
equations (5), (10), and (11) for the wavelength \ = 1.55 um, 
where water only weakly absorbs radiation, was investigated. 
Selected results are presented in Figures 3a—3c. It was 
assumed that the effective radius of droplets is equal to 6 
um and the parameter s = 6 in equation (3). The Mie 
calculations for this case yield: wo = 0.9935, g = 0.8214. It 
follows from Figures 3a and 3b that equation (30) is an 
accurate representation for tT > 10, 0 =0°, Vo = 60°. The error 
is less than 2% for this case, which is a relatively small error, 
compared to the uncertainty in forward cloud models applied 
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to the real atmosphere (e. g., vertically and horizontally 
inhomogeneous cloud fields). The constant error at tT > 30 
is mostly due to the error of the approximation for the 
reflection function ofa semi-infinite cloud, given by equation 
(32). The error increases with the optical thickness and 
reaches 5.5% at tT = 5 (see Figures 3a and 3b) in the case 
under study. 

[34] The error depends on the solar zenith angle. This is 
shown in Figure 3c. It is smaller than 5% for solar zenith 
angles, ranging from 12° to 75° at t = 10 and nadir 
observation. The error reaches 22% for the optical thickness 
T = 5 and the solar zenith angle 10°. It is, however, around 
5% at the solar angle of 60° for the same optical thickness 
(see Figure 3c). Thus, equations (30)—(32) can be applied to 
clouds with T < 10 with a great care. In particular, the solar 
angle should be in the range 50°—65° to have an accuracy 
lower than 7% at tT = 5. More complete studies of accuracy 
are given by Kokhanovsky and Rozanov [2003]. They also 
introduced a correction term in equation (30), which 
allowed to reduce errors of equation (30) well below 5% 
at t = 5 for all solar angles smaller 75°. However, we will 
not use the modified solution given by Kokhanovsky and 
Rozanov [2003] here and consider clouds with tT > 10. This 
allows to simplify the cloud retrieval algorithm presented 
here. The modified solution, given by Kokhanovsky and 
Rozanov [2003], can be easily incorporated in the cloud 
retrieval scheme, if clouds in the range of the optical 
thickness 5—10 are of interest. 

[35] The comparison of data for wavelengths 0.65 and 
1.55 um, presented in Figure 3a, show us that the limit of 
the semi-infinite cloud is reached more rapidly for infrared 
absorbing wavelengths, where water is a more absorbing 
substance than in the visible. As it will be shown in the next 
section, this result can be used in the estimation of the 
droplets size even if the optical thickness of clouds itself is 
not retrieved. 

[36] Interestingly, both curves in Figure 3a cross around 
the optical thickness 10. At optical thicknesses lower than 
10, the reflection function at the absorbing channel is 
higher. This is because of the differences in the phase 
function between the visible and infrared spectral regions. 

[37] Equation (30), which includes equation (6) as a 
special case, can be used to determine the optical thickness 
and effective radius from spectral reflection function meas- 
urements over extended cloud fields. 

[ss] The Lambertian surface reflection is accounted for 
by equation (12). The substitution of equation (30) into 
equation (12) yields: 


i At 


R(L, Hos) = Roo(H os) = }exp(—x — y) — T| T(P Ho): 
(33) 

where 
T (1, Ho) = tKo(})Ko(to) (34) 


is the transmission function of a cloud, r is the total 
reflectance of the cloud or the spherical albedo. The values 
of ¢ and Rẹ are given by equations (31) and (32) 
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respectively. Strictly speaking, equation (33) is only valid 
for channels, where the single scattering albedo wo is close 
to one. Note that for typical values of the cloud droplet 
effective radii 4—16 um the value of wo = 0.985—0.996 for 
the wavelength 1.55 um, which is indeed close to one. For 
the wavelengths 1.64 and 2.15 um, which are often used in 
cloud retrievals, the values of wo lay in the range 0.991— 
0.998 and 0.970—0.993 respectively for the same interval 
of the cloud effective radius change. Overall, the 
probability of photon absorption 8 = 1 — wọ is smaller 
than 0.035 for wavelengths \ < 2.3 um and aef < 16 um. 
In this case the parameter y < 1.1. 

[39] Equation (33) is a key formula for this study. Its 
accuracy (also with account for correction terms) is studied 
in detail elsewhere [Kokhanovsky and Rozanov, 2003]. 

[40] Note that equation (33) can be easily generalized to 
account for the bidirectional surface reflectance effects. We 
will not discuss this matter in any detail here, however. 


2.3. The Total Reflectance 


[41] The only parameter in equation (33), which is not 
represented by simple functions, is the total reflectance r. 
Let us find now the approximate solution for the value of r 
in equation (33) in terms of elementary functions. Clearly, 
the value of r 4 1 — t due to light absorption effects. 

[42] To start with, we remind that the total reflectance or 
the spherical albedo r is defined by the following equation 
[Sobolev, 1972]: 


n/2 
E= / do sin20oR(9o), 
0 (35) 


Qn n/2 
R(Vo) = I dọ J dð sin20R(90, 9, p) 
0 0 


Here, R(o) is the plane albedo [Sobolev, 1972]. For the 
case of idealized semiinfinite nonabsorbing clouds [Sobolev, 
1972] and as a result of the conservation of energy law, we 
have: 


20 n/2 n/2 
1 
x [af dòosin2 f dðsin2ðR? (o, 0, p) =1 (36) 
T 
0 0 0 
and also 
l Qn n/2 
x, | ef d§ sin 20R°. (90,9, p) = 1, (37) 
T 
0 0 


i.e., all photons injected into a cloudy medium are reflected 
back into outer space after an infinite travel time. Here, 
R? (ðo, Ù, p) is the reflection function of a semi-infinite 
nonabsorbing cloud. The reflection function Re (Do, D, p) 
of a cloudy medium only weakly depends on its 
microstructure (see Figure 1) and (by definition) it does 
not depend on either the optical thickness T = 0,,,/ or the 
single scattering albedo wo = Oscq/Oex Here og; is the 
extinction coefficient and osca is the scattering coefficient of 
a cloudy layer of the geometrical thickness L. 
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[43] It follows from equations (30) and (35) for light 
absorbing clouds: 


r = r% — texp(—x — y), (38) 
where 
n/2 
Too =| dvo sin 209 (Vo), 
(39) 
Qn n/2 
Ro(0o) = = I ao f dò sin 2R (ðo, , $). 
T 
0 0 
The normalization condition [van de Hulst, 1980] 
1 
2 fam Ko(u) = 1 (40) 


0 


was used to derive equation (38). Note that approximation 
(10) is consistent with equation (40). Now, we need to find 
the expression for the value of ra, given by equation (39). 
The constant rœ represents the total reflectance of a semi- 
infinite layer. According to the definition (36), 7, = 1 at wo = 
1. Equation (39) is not readily analytically integrated at 
arbitrary values of wo, however. So we make the approx- 
imate integration in equation (39). First of all, we note that it 
follows from equation (39) as wọ — 1 (see also the 
discussion in the previous section): 
Io = 1—y. (41) 
To derive equation (41) we have used formulae (18), (36), 
and (40). At larger values of y, we will use the same 
substitution as was used in the derivation of equation (28) 
from equation (25). Then we have for the integral (39): 
Yoo = exp(—y) (42) 
or (see equation (28)): ra = /'“. Note that it follows: a ~ 1 
and / % ræ. This is an interesting result by itself. The 
accuracy of equation (42) was studied in detail by 
Kokhanovsky [2002b]. It can be applied at least up to y = 1. 
[44] Combining equations (31), (38), and (42), we have 
for the total reflectance of a cloud layer: 


sinh(x + ay) exp(—y) — sinh(y) exp(—x — y) 
r= 
sinh(x + ay) 


. (43) 


The substitution of equation (43) into equation (33) enables 
us to obtain the reflection function of a cloud-underlying 
surface system in terms of the elementary functions, 
assuming that the reflection function of a semi-infinite 
nonabsorbing cloud is known (see equations (32) and (5)). 
This is the major result of this work. 

[45] Having r and ¢, one can also find the plane albedo 
[Zege et al., 1991; Kokhanovsky, 2001] R(uo) = Ralo) — 
T(uo)exp(—x — y), the diffuse transmission 7(\19) = tK (uo) 
and the absorptance A(j19) = 1 — R(uo) — T(uo). The diffuse 
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transmission is defined by an equation similar to equation 
(35) for a plane albedo, but with the transmission (34) and 
not the reflection function (33) as an integrand. Note that 
Ralo) = exp(—yKoluo)) is the plane albedo of a semi- 
infinite layer [Zege et al., 1991; Kokhanovsky, 2001]. 
Functions R(io) and T(uo) are often measured experimen- 
tally (see results presented by Los and Duynkerke [2000]). 
For visible channels, where light absorption can be 
neglected, we have: R(uo) = 1 — T(uo). Functions Kolpo) 
and T({19) decrease with the solar angle. This means that the 
plane albedo increases with the solar angle, which is 
opposite to the dependence of the reflection function at a 
nadir observation, given in Figure 1. Smaller values of the 
plane albedo for higher Sun positions were also observed 
experimentally [Los and Duynkerke, 2000]. 

[46] The total light absorptance in a cloud layer is given 
by: a=1 — r — t, where r is calculated by equation (43) and 
t by equation (31). As a result, we have the following 
analytical equation for the total light absorption inside a 
plane-parallel cloud with a finite optical thickness: 


_ sinh(x + ay) (1 


exp(—y)) — sinh(y)(1 — exp(—x — y)) 
sinh(x + ay) ; 


(44) 


where the diffuse illumination of a cloud layer is assumed. It 
follows as œ = 1 from equations (31), (43), (44): 


sinh(x) f sinh(y) 
pz — 
sinh(x + y)’ sinh(x + y)’ 
sinh(x + y) — sinh(x) — sinh( y) 
p= 
sinh(x + y) 


(45) 


? 


which yields the well-known formulae, presented elsewhere 
[e.g., Zege et al., 1991]. 


2.4. Similarity Parameters 


[47] We see that the global radiative characteristics of 
cloudy media are well described by only two parameters: 
x= 7,/3(1 —wo)(1—g) and y=4,/7744. We will call 
them similarity parameters. Let us discuss their relationships 
with microphysical properties of clouds, in particular with 
the sizes of droplets and their complex refractive index. 

[48] The physical sense of the parameter y was discussed 
above. It is equal to the logarithm of the inverse value of the 
spherical albedo of a semiinfinite weakly absorbing light 
scattering cloud layer (see equation (42)). The parameter x 
describes the attenuation of a light field in deep layers of 
semi-infinite weakly absorbing media. Namely, we have for 
the light intensity in deep layers: (u, T) ~ i(\)exp(—x), 
where the angular distribution of light field i(j1) does not 
depend on the optical depth t. Thus, we see that radiative 
characteristics of optically thick layers are determined by 
parameters, which govern light reflection and asymptotic 
regime for semiinfinite media. This is an important result of 
the radiative transfer theory, which has not yet been fully 
exploited for the solution of remote sensing problems. 

[49] The transformation of x and y yields another set of 
parameters: p = xy and q = x/y, where p = 4Tabs, q = 30%, 
Here Tabs = Oabsl is the optical thickness resulting from 
light absorption and T* = o,,(1 — g)L is the scaled optical 
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thickness, which is smaller than the cloud optical thickness 
by a factor (1 — g). The following parameters such as the 
absorption length Labs = l/Oabs and the reduced extinction 
length Lë = 1/[o..41 — g)] can also be introduced. There- 
fore, we see that values of r, t, and a for optically thick 
weakly absorbing media only depend on the ratios L/Labs, 
L/L*.,. The details of the angular distribution of light in a 
single scattering event are of a secondary importance for 
this assumed case. 

[so] Let us estimate the dependence of parameters p and q 
on the effective radius aef of droplets. For this we will 
assume that droplets are large (kagy — oo) and weakly 
absorbing (Ka. — 0), which is appropriate for water 
droplets in the visible and near-infrared regions of the 
electromagnetic spectrum. Here k = 21/) is the wave 
number and k = 4nx/ is the absorption coefficient of 
water (x is the imaginary part of the refractive index of 
droplets and A is the wavelength of the incident light). To 
the first coarse approximation (see Appendix B), extinction 
and absorption coefficients are given by: 


1.5C, 
—, Cabs = 1.25KC,,1 — g = 0.12 


aef 


Dai = (46) 
where C, is the volumetric concentration of droplets. The 
value of C, is the fraction of a unit volume occupied by 
droplets in a cloud. It is a small figure. For example, assuming 
that C, ~ 1077 we obtain: o ex œ 0.0257! at ar 6 um, which 
is inside the variability range for the value of dew in natural 
clouds [Zege et al., 1991; Los and Duynkerke, 2000]. 

[sı] The probability of photon absorption B = Oabs/Oext 
can also be obtained from equation (46): 8 = 2 Kder- For p 
and q we have to a first approximation, therefore: 

(47) 


p=Skw/p, q=0.14w/pae, 


Here w = C,Lp is the liquid water path and p = 1000 kg/m? 
is the density of water. 

[52] This means that the parameter p does not depend on 
the size of droplets at all in the framework of the approx- 
imation considered. The value of p is determined only by the 
product of the absorption coefficient of water and the liquid 
water path of a cloud. The parameter g, on the other hand, 
does not depend on the refractive index of droplets. It is 
determined by the ratio of the liquid water path w = C,Lp to 
the effective size of droplets aep This means that both the 
imaginary part of the refractive index of droplets x and the 
values ay, w can be in principle determined from the measure- 
ments of the backscattered radiation from clouds. 

[53] In the next section it will be assumed that the spectral 
dependence k(\) is known. This assumption allows us to 
find the effective radius of droplets and the liquid water path 
of cloudy media from satellite measurements of the reflec- 
tion functions of clouds. Note, however, that the uncertainty 
related to imperfect knowledge of the temperature-depend- 
ent absorption coefficient k(\) for real cloud systems could 
be quite large [Asano et al., 2000]. 

[54] Using relationships between (x, y) and (p, q) one 
obtains approximately: 


4 
r= S y = 6 Kae, (48) 


z ap 
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where the value of y is determined by the effective size of 
droplets and x ~ w/,/a. Although the dependencies 
presented in equation (48) are only approximations of a real 
situation in cloudy media, they illustrate the origin of the 
parameters used to describe the local optical characteristics 
of clouds quite well. It should be pointed out that the 
reduced optical thickness T* = (1 — g) (see equation (7)) is 
given by the ratio 4x/3y or T* © 8w/45paey: 


3. The Inverse Problem 
3.1. The Cloud Optical Thickness Determination 
[ss] Now all equations are prepared to proceed to the 


solution of the inverse problem. In particular, parameters aef 


and w can be found as roots of following functions: 


F; (der, w) = Ri(ut, Ho, 9) — Rooi (H Ho, 9) 
te = Ait; (tj + re) 
1- Airi 


Ko(uo)Ko(H), (49) 
where i= 1, 2, 3, ...N corresponds to a given measurement 
channel and R; is the measured reflected function at that 
channel. N is the total number of channels used in a retrieval 
procedure. For reasons of simplicity and clarity, we assume 
that all possible atmospheric influences (like gaseous and 
aerosol absorption and scattering) are subtracted from 
measured reflection functions and R; is the so-called refined 
reflection function. For the same reason, we avoid the 
integration of equation (49) with respect to the spectral 
response function of the instrument and the spectral change in 
the solar constant for a given instrument response function. 

[s6] In this paper we will only consider one-channel (N = 
1) and two-channel (N = 2) algorithms. The generalization to 
the case of an arbitrary number of channels is trivial. The two- 
channel algorithm constitutes a minimum set of experimental 
values to determine two unknown parameters (aep and w) as 
simultaneous roots of functions F (aep w) and F(a, w). This 
does not mean, however, that the one-channel algorithm does 
not provide us with any useful information. To be more 
explicit, let us suppose that the wavelength, of which the 
measurement is performed, belongs to the visible range of the 
electromagnetic spectrum. In this case we have with a high 
accuracy (see equation (47)): p = 0 and we are left with the 
necessity of finding only one parameter q = x/y, which is the 
root of the function F'\(qg). Then the value of T* = tq can be 
easily found. 

[57] The discussion of the two-channel algorithm will be 
postponed to the next subsection. Let us consider the one- 
channel algorithm (at wọ = 1) in more detail. For this we 
need to return to equation (14), which gives: 


(1—A)A 


“TA ay Po 


where the function A is introduced and given by 


R? (u; Hos 9) — RCs Ho» P, T) 

A = A(t, 9) = 22 i 51 
(epote) Ko(u)Kalt) ae 

The analytical results for functions R? (p, Ho, 9) and Ko() 
have been presented above (see equations (5) and (10) and 
Appendix A). Thus, the global transmittance ¢, and, 
correspondingly, the total reflectance r = 1 — £ of extended 
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cloud fields can be obtained from equations (50) and (51), 
and a knowledge of the surface albedo A, angles Jo, 0, p and 
the measured value R(t, uo, Y, T). For such a retrieval, one 
does not need to know the optical thickness of clouds and the 
average size of droplets. This is of significance for climate 
studies, where the globally and temporally averaged value of 
the cloud reflectance r is an important parameter. Note that 
for measurements over oceans (and over land for sufficiently 
small values of A (e.g., 400 nm)), the value of A is small and 
can be neglected in a first coarse approximation. 

[ss] Let us consider equation (50) in more detail now. 
First of all, we have at A = 0: t= A and, second, t= 0 and r 
= 1 at A= 1. This shows that all photons entering optically 
thick clouds over surfaces with A = 1 survive and return 
back to outer space. They yield no information about actual 
cloud thickness. This explains why the retrieval of cloud 
parameters over bright surfaces (e.g., snow and ice) cannot 
be performed in visible. This issue, however, will be 
discussed in the next section in more detail. 

[s9] The information on the global transmittance ¢ can be 
used to find the scaled optical thickness T* = T(1 — g). It 
follows from equation (7): 


(52) 


where f is given by equation (50). Again the value of T* can 
be obtained without knowledge of the size of droplets and 
the actual optical thickness of clouds. 

[60] Equation (52) can be used for the retrieval of t* from 
the measurement R(p, po, p, T) at a single wavelength. The 
functions R? (u, po, ~) and Kolpo) in equation (51) and the 
parameter & in equation (52) are defined by equations (5), 
(10), and (11), respectively. 

[61] One can apply equation (50) for the derivation of the 
optical thickness t as well if the value of g is known. 
However, the asymmetry parameter g only slightly depends 
on the size of droplets for nonabsorbing channels [Kokha- 
novsky, 2001]: 

= £0 - ———; 53 

& = 80 aa (53) 
where go ~ 0.884, Y ~ 0.5, k = 4. Often, the dependence 
g(a) is neglected and it is assumed that aep = 10 um 
[Rossow and Schiffer, 1999]. In this case it follows from 
equation (53): g = 0.86 at \ = 0.65 um and a,¢= 10 pm. This 
value of g can be used for a crude estimation of the optical 
thickness from the value of T*. 

[62] Clearly, errors can be introduced, if one assumes the 
fixed a priori defined value of g. Indeed the value of aep is 
usually in the range 4—20 um [Han et al., 1994]. Then it 
follows from equation (53) at \ = 0.65 um that g = 0.84— 
0.87 at ae = 4—20 um. Also we have: tT = ut*, where v = 
(1 — g) | © 6.3-7.6 and TE [9.4, 11.5] at T* = 1.5, 
depending on the value of g used. The assumption that aef = 
10 um yields: g = 0.86 and v = 7.2, t = 10.7. This leads to a 
relative error of 7—14% in the retrieved optical thickness, 
depending on the actual value of a.» This uncertainty in the 
optical thickness determination can be removed if measure- 
ments in the near-infrared region of the electromagnetic 
spectrum are performed, enabling the size of droplets and, 
therefore, the parameter g (see equation (53)) to be estimated. 


AAC 4-12 


3.2. The Two-Channel Inversion Algorithm 


[63] As it was specified above for the correct estimation 
of the optical thickness of clouds from space we need to 
know the effective radius of droplets. The size of droplets 
can be found if the reflection function in the near-infrared 
region of the electromagnetic spectrum is measured simul- 
taneously with measurements in the visible range. This is 
due to the fact that the reflection function in the infrared 
strongly depends on the probability 8 of photon absorption 
by local volume of a cloud medium. This probability is 
proportional to the effective radius of droplets, as it was 
discussed above (8 ~ kaef) (see section 2.4). 

[64] The influence of absorption and scattering of light by 
molecules and aerosol particles on the measured value R(p, 
Ho, Y, T) is neglected in the retrieval algorithm presented 
here. However, corrections can be easily taken into account 
if needed [Bucholtz, 1995; Wang and King, 1997; Goloub et 
al., 2000]. The influence of the surface reflection on the 
cloud reflection function, assuming that the surface is 
Lambertian with albedo A, is taken into account. Then it 
results (see equations (14) and (33)): 

Ay (ag, w) = wo, -iC Al 

1 —Ai[1—t (ae, w 


y] Ko()Ko(tto), (54) 


Ro (aep, w) = RY, exp(—y (aef) (1 — cy (aes) )u) 


h (der, w) Ay 
1 — Aor (ae, w) 


an exp(—x(aey, w) y(ae)) 


` t (aef, w) Ko (1) Ko (Ho). (55) 


[6s] The subscripts “1” and “2” refer to wavelengths , 
and ^z in visible and near-infrared channels, respectively. 
The values of 4, and A> give us the surface albedos in 
visible and near infrared, respectively. We assume that both 
A, and A» are known. R; and R, are measured reflected 
functions at the visible and near-infrared channels, respec- 
tively. Any weak spectral dependence of the function RO) 
due to the third term in the nominator of equation (5) is 
considered negligible because of small values of this term 
for most geometries. The explicit dependence of functions 
involved on the parameters ae and w we are going to 
retrieve is introduced in brackets. The liquid water path w 
is preferred to the optical thickness in the retrieval proce- 
dure due to the independence of w of the wavelength. The 
optical thickness is uniquely defined if a,-and w are known 
(see Appendix B). 

[66] Equations (54) and (55) form a nonlinear system of 
two algebraic equations, having two unknowns (aef and w). 
Standard methods and programs are available for the 
solution of such problems. However, taking into account 
the specific form of equation (54), the exclusion method is 
used. This allows a single transcendent equation with one 
unknown (ae) to be formulated. It follows from equation 
(52) that 


A(t" (ag) - 9) 
3(1— gı (aer)) ` 


[67] On the other hand, the optical thicknesses Tı is 
related to the liquid water path w and the effective radius 


(56) 


TS 


KOKHANOVSKY ET AL.: SEMIANALYTICAL CLOUD RETRIEVAL ALGORITHM 


of droplets by the following equations (see Appendix B): 


q= Woext (An; aef); (57) 
where Gey = o3/p, p is the density of water and o%, is given 
by equation (B3). We have from equations (56) and (57): 


w- p(t! (ae) = a) 
30% r,ag) (1 — g1 (ag) ) 


where the dependence of functions on the value of the 
effective radius is explicitly noted. The value of ¢; is obtained 
from equation (54) (see also equations (50) and (51)). Thus, 
one unknown parameter (w) is expressed in terms of another 
parameter (aef). 

[6s] The substitution of equation (58) into equation (55) 
yields a single transcendent equation for the determination 
of the effective radius of droplets in a cloud. The equation 
obtained can be easily solved numerically. In particular, we 
use the standard subroutine zbrent described by Press et al. 
[1992]. This subroutine relies upon the Brent’s method of 
finding roots of the function [Brent, 1973]. The obtained 
value of ay is used to calculate the optical thickness via 
equation (56). The liquid water path is determined from 
equation (58). This concludes the two-channel algorithm 
description for the determination of cloud parameters. 

[69] Note that the subroutine zbrent performs the root 
search with a high speed. In practical terms, the speed is 
close to the speed of calculating an expression, consisting of 
ordinary elementary functions. The results are obtained 
immediately in real time. In this sense, we can drop the 
word “semi” from the title of our paper and treat our 
algorithm as a truly analytical inversion procedure, which 
is a surprising result, taking into account well-developed 
multiple light scattering in a cloud and large size of droplets. 
On the other hand, namely these two conditions allowed us 
to develop the algorithm. 

[70] The algorithm proposed is extremely simple, robust, 
and fast. Moreover, the insight in the factors involved and 
equations presented can be of help for the general cloud 
optical studies, not necessarily related to the inverse prob- 
lem. It could be also of interest to compare results of our 
retrieval procedure with other methods (and simultaneously 
with in situ measurements) for real cloud fields. We expect 
that, owing the complexity of cloud fields and all uncer- 
tainties involved [Pincus et al., 1995], the method presented 
here will have in the end an accuracy similar to that 
obtained using the exact radiative transfer equation. This 
is mostly due to the fact that the notion “exact” is often (if 
not always) not relevant in studies of complex time and 
space-dependent fields of geophysical parameters. 

[71] We would like to make one more remark. In some 
cases it is highly desirable to minimize the influence of 
generally unknown and potentially significant surface 
reflection. This is achieved by making both measurements 
in the infrared region of the electromagnetic spectrum 
[Platnick et al., 2001]. Then equation (54) becomes: 


(58) 


Ait 


R, =R? 1 
1 = Ry exp(—yi( 1 An 


cy )u) exp(—x1 — yı) 


- t1Ko()Ko(Ho). (59) 
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where values of t, 71, x1, yı have the same meaning as in 
equation (55), but for another wavelength. The only 
difference is that the wavelength ^; is now in the near- 
infrared range of the electromagnetic spectrum. 

[72] The retrieval procedure is based on the solution of 
two nonlinear equations (55) and (59) in this case. This can 
be easily performed numerically. It should be stressed, 
however, that the sensitivity of the reflection function in 
the infrared to the liquid water path is low. This is due to the 
fact that the function 


Q(p, Ho; P, T) = Ræ(p, Ho, P) — Îi (p, Ho, P, T) (60) 


approaches zero with the increase of T more rapidly in the 
infrared than it does in visible (see Figure 3a). This also 
means that the second term in equation (55) can be omitted 
in this case. This allows us to find the value of the cloud 
droplet effective radius for optically thick clouds with 
unknown values of the cloud optical thickness. 


3.3. The Error Propagation Study 

[73] The influence of various possible errors and uncer- 
tainties on the retrieved values of aep and T needs to be 
investigated. This was done by many authors [King, 1987; 
Rossow et al., 1989; Pincus et al., 1995; Asano et al., 2000]. 
However, for the sake of completeness, we need to revisit 
this issue. 

[74] Let us consider first of all the one-channel algorithm, 
where all calculations can be done analytically [King, 
1987]. In particular, we have the following result for the 
optical thickness as a function of the measured reflection 


function R: 
1 fu A 
j ea. eu 
where 
p= tieg (62) 
and 
R 
v= a (63) 


We have used equations (56) and (50), omitting index “1” 
for the sake of simplicity. The viewing function u is defined 
by equation (29), A is the surface albedo, and a is given by 
equation (11). It follows from equation (61) that the 
uncertainty in T can be presented by the following equation: 


OT OT OT OT 
At= 7,44 | ay a” | au A" | ap Ab (64) 
where 
OT z 1 OT u ðr 1 OF T (65) 


OA b(1 — A)?’ Ov by?’ Ou bv’ Ob b` 
Similar equations were derived earlier by King [1987]. The 
analysis of equations (64) and (65) allows to make the 
following conclusions: 

1. The error of the determination of the optical thickness 
increases with the value of albedo as (1 — A)”. In 
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particular, or =œ as A — |. This implies that it is 
impossible to find the optical thickness of any optically 
thick nonabsorbing cloud at A = 1. The physical grounds for 
such an effect are quite simple and will be discussed now. 
All photons transmitted by the cloud are reflected back to 
the cloud by the reflection from the underlying surface and 
then they return back to outer space. Thus, the value of the 
reflection function is that of a semiinfinite medium, which 
by definition does not depend on the optical thickness. The 
same result follows from equation (14) at A = 1. This makes 
it extremely difficult to retrieve the cloud optical thickness 
over white surfaces (e.g., snow, ice, and underlying cloud 
layers with high light reflectance) in the visible [Platnick et 
al., 2001]. To avoid this problem one should make 
measurements in infrared or even microwave regions of 
the electromagnetic spectrum, where the albedo of natural 
underlying surfaces is low and atmospheric emission is 
high. It should be pointed out, however, that the value of the 
reflection function as a function of the optical thickness in 
the infrared region of the electromagnetic spectrum reaches 
a saturation value quite rapidly due to the high values of 
light absorption (see Figure 3a). This also makes it difficult 
to retrieve the optical thickness of thick clouds from 
measurements in infrared channels. 

2. We have as T — œo: R — R» and v — 0. Note that 
the uncertainty in the retrieved value of the optical 
thickness due to the uncertainty in the value of v increases 
as v ~ as v — 0. This makes it impossible (see equation 
(65) for or) to derive the optical thickness of extremely 
thick clouds with T > Tmax- The upper limit of the optical 
thickness Tmax depends on many factors, including the 
surface albedo A, the maximal relative error of measure- 
ments 6, the observation geometry, the solar zenith angle 
and the value of the asymmetry parameter g. Let us define 
the value of the upper theoretical boundary of the optical 
thickness Tmax» which can be retrieved in principle from the 
reflection function measurements in the optical range by the 
following equation: {(Tmax) = ô, where the global transmit- 
tance f is given by equation (7). Clearly, the value of T 
cannot be obtained for T > Tmax from reflection measure- 
ments in visible. Choosing the value of g = 0.85 and the 
error of measurements ô = 0.05, we obtain that Tmax = 168. 
The actual optical thickness of clouds is usually within the 
interval of 5—120 with the most frequently occurring value 
of around 30 [Trishchenko et al., 2001]. So one can assume 
that the optical method of the cloud optical thickness 
determination does not have limits for large T in practical 
terms. However, this is not true [Pincus et al., 1995]. In 
particular, our estimation of Tmax was performed for black 
underlying surfaces. Actually, the maximal optical thickness 
Tmax depends on the surface albedo. It is smaller for larger 
surface albedos. It is smaller at smaller solar zenith angles 
(at nadir observation) due to larger contribution of multiple 
scattering in this case [Pincus et al., 1995]. Actually, our 
estimation of Tmax should be substituted by Tmax — AT 


where [Zege et al., 1991] At = a for a nonzero 


surface albedo. We see that the actual Sale of the maximal 
optical thickness is negative at A close to 1. So it cannot be 
retrieved as it was discussed before. Values of A could reach 
0.97 for snow surfaces. Then we have, assuming that g = 
0.85: At = 116. This means that the actual maximal cloud 


optical thickness, which can be retrieved over snow areas in 
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visible is around 50. Actually, due to other uncertainties 
[Pincus et al., 1995], it is even smaller. 

3. Clearly, one can use infrared measurements in this 
case. The albedo of the snow surface is quite low in the 
infrared. However, another problem arises: light absorption 
by water cloud itself, which increases the cloud reflection 
function saturation rate. So one should use even larger 
wavelengths to avoid problems, which occur in the optical 
band (e.g., microwave cloud sounding from satellites). 

4. The uncertainty in T due to the uncertainty in the value 
of u is proportional to v™! (see equation (65)). Therefore, it 
increases as v — 0. This makes it important to have an 
accurate approximation for the viewing function u (see 
equation (29)) for the case of thick clouds. 

5. The uncertainty in T due to the uncertainty in the value 
of b is proportional to t/b (see equation (65)). Thus, it is 
more important for larger values of the reduced optical 
length t/b than for smaller ones. Larger values of g result in 
larger uncertainties in the optical thickness determination. 

[75] Let us consider now the uncertainties in the retrieved 
parameters in the two-channel algorithm. Unfortunately, it is 
more difficult to do analytically. So we study this case using 
the numerical approach. The reflection function of a water 
cloud having a.¢= 10 um, T = 10, Vo = 49°, 0 = 7°, p= 0° 
and )\; = 0.65 um, 2 = 1.55 um is first calculated with the 
doubling-adding radiative transfer code. Then we introduce 
the error of measurements 6 in both channels and retrieve 
values of a. T, using the algorithm, described in the previous 
section. Results are presented in Figure 4. As expected, the 
error of retrieval A increases with the value of ô. The value of 
A is not zero at ô = 0 due to the approximations made in the 
formulae used in the inversion scheme. The error of retrieval 
is less than 10% for A = 5%, which is not untypical for 
satellite radiometric measurements. Note that the observa- 
tion angle in this example is not exactly zero. However, 
the accuracy of retrieval is quite high. It suggests that 
equation (5) with bı = 1.48, b = 7.76 can be used at other 
than nadir observation angles, providing that cos Jo ~ 1. 
The modification of equation (5) to account for arbitrary 
observation geometries will be considered in a separate 
publication. 

[76] The second example is presented in Figure 5. For this 
case the value of ô is assumed to be equal to 5% and the 
influence of the optical thickness on the retrieval accuracy is 
investigated. One can see that the retrieval of t for the case 
of larger optical thicknesses is more difficult. This is due to 
the reflection function being nearly independent of the 
thickness of a cloud at large t (the limit of a semi-infinite 
medium) as discussed above. 

[77] The error of the effective size retrieval initially 
increases with the value of t. However, it approaches a finite 
value of approximately 20% independently of the optical 
thickness tT. This result is also easily understood. Indeed, 
equation (55) transforms into equation (32) as T — oo. The 
latter does not depend on the optical thickness. The error in 
this case can be retrieved analytically from equation (32). 
This is due to the fact that the two-channel algorithm is 
reduced to a one-channel algorithm (but for the effective 
radius of droplets determination). Namely, we have: 


1 es 
=== [In— 
y RO 


u o0 


(66) 
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where we neglected the value of c for the sake of simplicity 
(see equation (32)). On the other hand, it follows from 
equation (48): 


y= Eyra, (67) 


where = = 6. Thus, it follows from equations (66) and (67): 


1 in? Ree 


r = — 68 
Aef Kn? n Ro’ ( ) 


where 1 = Zu. Therefore, we obtain for the uncertainty in 
def: 


Oa, Oa 


Oder ef ef 
Aas = — An + — An AQ, 69 
aga a Aha, Ant Bae (69) 
where O = R°./R,, and 
Oder ag Ode 2p Oas 2Gey (70) 
On k’ On n ` Q Old 


From equations (69) and (70), the following can be 
summarized: 

1. The relative error of the retrieval of the effective radius 
is equal (with opposite sign) to the error with which the 
absorption coefficient of water is known (at An = AQ = 0). 
The value of Ak/k is of the order of several percent, 
depending on the a priori unknown temperature of droplets 
[Kuo et al., 1993; Asano et al., 2000]. Note that the value of 
k can be also modified by the presence of soot or aerosol 
particles trapped inside a droplet. This is a typical situation, 
e.g., over large industrial urban areas and cities. Obviously, 
the concentration of impurities also plays a role. This might 
well explain a discrepancy (constant shift of 5%) between 
satellite retrieved data for the effective radius and that 
obtained from in situ measurements [Nakajima et al., 1991]. 
In situ measurements were performed using the laser 
diffractometer [Nakajima et al., 1991] and the results, 
therefore, are largely insensitive to the value of the 
absorption coefficient. This is due to the fact that the 
information on the size of particles is obtained from 
the analysis of the diffraction pattern. This pattern almost 
entirely depends on the geometrical characteristics of 
droplets (primarily, on their diameter or droplet size 
distribution in the measured volume). This is, of course, 
not the case for the reflection function, which depends on 
the product Kaer (see equation (48)). 

2. The relative error of the retrieval of the effective 
radius of droplets is two times larger than (but with 
opposite sign) the error of the function ņ = Eu (at Ak = 
AQ = 0). The error of the function ņn = Ew results 
primarily from error in the approximate equation for the 
function R&(\1, uo, 9), which is approximately of 2%. This 
alone yields a uncertainty in the retrieved value of the 
effective radius of 4%. 

3. It follows for the relative error of the effective radius 
determination at Ak = Ay = 0: 


Aag _ p AQ 


Aef Q’ 


(71) 
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Figure 4. The dependence of the error of retrieval of a,-and T on the error of measurements A at Jo = 49°, 
0 = 7°, p = 0°, A; = 0.65 um, ` = 1.55 jum. It was assumed that a,-= 10 um and T = 10 at A = 0.65 pm. 
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Figure 5. The dependence of the error of retrieval of a.yand Tt on the optical thickness T at a.-= 6 um, 
Ùo = 49°, 0 = 7°, 0 = 0°, A; = 0.65 um, Az = 1.55 um, and ô = 5%. 
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where P= 2 In” 'Q is the error amplification coefficient. This 
coefficient is equal to infinity at Q = 1 as one can expect. It 
follows from our calculations that Q = 1.5 at a.¢= 6 um, \ = 
1.55 um, Ùo = 49°, 0 = 7° and ọ = 0°. This means that P ~ 4 
for this case and, correspondingly (see equation (71)) 
A2 L 5% at 92 = 5%, which is close to the asymptotical 
value, determined numerically (see Figure 5 for the value of 
ĉa ag T —+ co). This confirms our calculations. 

“T78] In our discussion we completely omitted other sour- 
ces of uncertainties like the imperfect knowledge of the 
atmospheric state, the in principle unknown cloud thermo- 
dynamic state, and the influence of possible subpixel cloud 
variability. Excellent reviews of these and other relevant 
issues, however, are available elsewhere [Rossow et al., 
1989; Kobayashi, 1993; Cahalan et al., 1994; Pincus et al., 
1995; Wang and King, 1997; Asano et al., 2000]. 


4. Conclusions 


[79] Modern satellite cloud retrieval algorithms use either 
the lookup table approach or the classical asymptotic theory. 
The classical asymptotical theory is valid only for optically 
thick light scattering layers (tT > 10). Clearly, both 
approaches converge as T — oo. In particular, King 
[1987] has shown that the difference between forward 
model results produced by two methods is below 1% at 
T > 9. This difference is completely negligible for all 
practical purposes. Note that most of the radiative transfer 
solvers also have numerical errors, which may reach 1% in 
some particular cases. 

[so] Thus, there is no need to apply the exact radiative 
transfer equation and, therefore, conventional lookup table 
approaches to the cloud satellite retrievals for cases of 
optically thick cloudy media. For optically thin clouds, 
however, the application of the lookup table approach is 
not avoidable. 

[sı] The classical asymptotic theory offers a great sim- 
plification of the cloud inverse problem solution and should 
be used in the analysis of all scenes where optically thick 
clouds are present. It should be stressed, however, that 
asymptotic formulae (see equation (15)) do not provide 
the analytical solution of the radiative transfer equation in 
terms of elementary functions. Their application is not 
simple. In particular, one has to solve a number of integral 
equations to find auxiliary functions and parameters 
involved in the main result of the asymptotic theory, given 
by equation (15) [Nakajima and King, 1992]. 

[s2] Asymptotic equations are valid at the arbitrary value 
of the single scattering albedo. The question arises whether 
it is possible to simplify these equations to account for a 
weak absorption of light in real terrestrial clouds (at least up 
to A = 2.3 um). The answer to this question is positive. To 
show this was the main idea of this paper. We also show the 
possibility to apply the developed scheme to the cloud 
inverse problem solution. Of course, obviously, this is not 
the only area of applicability of the equations obtained. We 
are currently working on generalizations of these equations 
for the case of inhomogeneous media and non-Lambertian 
surface reflectances. Other applications (e.g., to ice, snow, 
and foam optics) are also possible. Detailed studies of the 
accuracy of equations derived and their generalizations for 
smaller t are given elsewhere [Kokhanovsky and Rozanov, 
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2003]. In particular, the error is less than 3% at tT > 8 for 
cases presented in Figure 3a. 

[ss] Our method has a generally lower accuracy than a 
conventional asymptotic theory method as far as a forward 
modeling is concerned. However, it does not automatically 
mean that it will produce larger retrieval errors dealing with 
measurements performed in the terrestrial atmosphere. 
Indeed, errors of the modified asymptotic theory appear to 
be smaller than measurement uncertainties, not to mention a 
fundamental problem of partially cloud covered pixels or 
pixels with inhomogeneous cloud fields. This problem will 
be of a special importance, for a Scanning Imaging Absorp- 
tion Spectrometer for Atmospheric Chartography (SCIA- 
MACHY), launched on 1 March 2002 on board of the 
ENVIronmental SATellite (ENVISAT) [Bovensmann et al., 
1999]. The algorithm described here is specifically designed 
for this instrument, which has an extremely wide field of view 
(approximately, 30 x 60 km, depending on the wavelength). 
Obviously, in this case uncertainties other than differences 
between modified and conventional forms of asymptotic 
solutions will dominate the retrieval accuracy. We plan to 
discuss these issues, however, in our next publication. 


Appendix A: The Phase Function 


[s4] Phase functions of water clouds at A = 0.65 um are 
presented in Figure Al. They have been obtained assuming 
the gamma particle size distribution 


= Dao e244 : 


f(a) (Al) 
where D = const and the effective radius aep was changed 
from 4 to 20 pm. It is a complex task to model curves in 
Figure Al with simple analytical equations. 

[ss] It should be stressed, however, that for most of cases 
the third term in the nominator of equation (5) (p(8)) is 
much smaller than the first and second terms. Thus, even the 
comparatively large error in the phase function p®) will not 
influence the reflection function Re. (0, Ùo, p) in equation 
(5) very much. This follows from Figure 1 as well. Thus, we 
will neglect the dependence of the phase function on the 
size of droplets and use the following approximation, 
obtained by the parametrization of computations with the 
exact Mie theory: 


(A2) 


5 
aay: 5 bje 88-8) 
i=l 


where @ is the scattering angle in radians, € = 17.7, o = 3.9 
and constants b; B; 6; are presented in Table Al. The 
accuracy of this approximation can be estimated from 
Figure Al. 


Appendix B: 
Characteristics 


Integral Light Scattering 


[s6] Integral light scattering and absorption characteristics 
of water clouds can be found with following equations 
[Kokhanovsky, 2001]: 

(B1) 


— raž — až 
Set = Ort Cv, Cabs = OäbsCy, 
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Figure Al. Phase functions of water clouds at different values of aep obtained with the Mie theory and 


approximation (A2). 


1 — g = 0.12 + 0.5 (kag) 70.1 5Kaey, (B2) 
where k = 27/A, A is the wavelength, a,;is the effective size 
of droplets, k = 4nx/, x is the imaginary part of the 
refractive index of water C, = N(V), is the volumetric 
concentration of droplets, N is number of droplets in a unit 
volume of a cloud, (V) is the average volume of a droplet in 
a unit volume of a cloudy medium, g is the asymmetry 
parameter, o,,, and Oaps are extinction and absorption 
coefficients. The values of o%, and o%*,, are given by 
[Kokhanovsky, 2001]: 


It follows from equations (B1) and (B3) for the ratio of 
extinction coefficients at two wavelengths: 


2/3 
LIS 
pstaa (B5) 
lil+¢ 
Table Al. Parameters b;, B; and 0; 
i bi B; 0; 
1 1744.0 1200.0 0.0 
2 0.17 75.0 2.5 
3 0.30 4826.0 T 
4 0.20 50.0 T 
5 0.15 1.0 T 


where, ¢ = an ,j = 1,2. Note that we have for large values 


of G — œ: ® — 1. 
[87] The optical thickness is given by: 


T = Onl, (B6) 
where L is the geometrical thickness of a cloud. It follows 
from equations (B6), (B1), and (B3): 
T = GeyW, (B7) 
where w = pC,L is the liquid water path, p is the density of 
water Gex = o%/p. Note that values of o,,, and Ož 
numerically coincide if afis expressed in jum and p in g/m’. 
Then w in equation (B7) is given in g/m°. The accuracy of 
equations (B1)—(B4) has been studied by Kokhanovsky and 
Zege [1995, 1997] and Kokhanovsky [2001]. It is better than 
5-8% at \ < 2.2 um. 
[ss] If higher accuracy is needed one can use the direct 
parametrization of the Mie calculations. We obtain with the 


Table B1. Parameters c; for Different Wavelengths 


A (um) Co c C2 C3 C4 
0.6457 0.1121 0.5118 0.8997 0.0 0.0 
0.8590 0.1115 0.4513 1.2719 0.0 0.0 
1.239 0.1095 0.4198 1.5796 0.0 0.0 
1.549 0.0608 2.465 —32.98 248.94 —636 

A (um) do dı dz ds d4 
1.239 1.779 —0.0068 4.94E—5 —1.425E-7 0 
1.549 1.671 0.0025 —2.365E—4 2.861E—6 —1.05E-8 
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accuracy better than 1% at selected wavelengths, which are 
almost free of gaseous absorption: 


4 


i+g=) c,(kag) ””, (B8) 
n=0 
ATX E n 
Oabs = Ea Cy 5 dn (kaes) , (B9) 


n=0 


where constants c,, d,, are presented in Table B1 at \ = 
0.6457, 0.8590, 1.239 and 1.549 um. There is no need to 
modify equation (B3) due to the high accuracy of this 
equation in the spectral range under consideration. 

[s9] At larger wavelengths it is easier to make paramet- 
rizations of parameters y and z = x/w instead of d ext O abs and 
1 — g. In particular, it follows at \ = 2.3 um: 


y = 0.24206 + 0.02524xer — 1.4197 - 10-*x2,, (B10) 
z = 0.02008 — 0.54709x,,"° + 7.53062x,  — 15.1524x,?. 
(B11) 


Note that the value of y is dimensionless. The dimension of 
z coincides with that of the liquid water path w. 
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